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Foreword 


This paper constitutes the- final technical report 
on a study of Surface: Properties of Ocean Fronts for the 
National Aeronautics and Space Administration under Contract 
No. NAS 5-22370. The basic report provides background infor- 
mation on oceanic fronts and describes the results of the 
several models which were developed to study the dynamics 
of oceanic fronts and their effects on various surface proper- 
ties. The details of the four numerical models used in 
these studies are given in separate Appendices; which contain 
all of the physical equations, program documentation and 
running instructions for the models. 

The following. scientists participated as consultants and 
scientific collaborators in this -project : Professors D.F. 

Liepper and R.J. Renard from the U.S. Naval Postgraduate 
School and Dr. Taivo Laevastu and Messrs. S. Larson and K. 

Rabe from the U.S. Navy' s Environmental Prediction Research 
Facility. 
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1.0 DEFINITION OF OCEAN FRONTS 

Although oceanic fronts, especially major ones, have 
been referred to in scientific literature at least as long 
as their atmospheric counterparts, their definition is not 
as precise and definite as for atmospheric fronts. This is 
partly due to the fact that the nature of oceanic fronts 
can be different in different locations and is affected by 
a great number of factors, both permanent and temporal ones. 
There has been .recently a renewed interest in the study of 
oceanic surface fronts, partly due to their importance in 
naval, fisheries, and navigation problems and partly due 
to the fact that most oceanic fronts can be observed from 
satellites, either as temperature or color or sea surface 
roughness discontinuities. 

The best definition of oceanic fronts is that of LaFond 
(1960), "The oceanic front is the leading edge of a border 
separating unlike water masses. Fronts can occur not only 
between water masses of different salinity but also between 
those differing in other properties, such as temperature". 
Before describing the characteristics of the front, we have 
to define "surface water types". The term "water, masses " 
refers to deep water masses (below the thermocline) where 
temperature and salinity can be considered as conservative 
properties and which are defined on temperature-salinity-: 
diagrams. The term "water type" has been applied in the 


1-1 



past to a singular point in the temperature-salinity 
diagram. Temperature and salinity are not as conservative 
in surface layers; however, surface waters can be charac- 
terized by other properties, such as biological (e.g., 
plankton) content, color, turbidity, and range of annual 
change of temperature and salinity (Schott, 1935; Laevas tu , 
1963). The term "surface water type" has been introduced 
to designate surface waters- with . definite but .relatively . 
uniform characteristics. 

Thus the oceanic .fronts at the-surf ace are boundaries 
of unlike surface water types; they form limits (boundaries) 
of natural regions of the oceans (Laevas tu and LaFond, 1970) . 
The fronts can be major ones, easily recognizable from re- 
latively pronounced changes of the properties (e.g. , tempera- 
ture) or they can be relatively broad transition regions with 
gradual: ..changes . 

LaFond and LaFond (1971) found that the California 
Front, which can be considered as an "average" front, was 
characterized by strong horizontal temperature (and salinity) 
gradients, large vertical changes in isotherm depths, tem- 
perature inversions , and a large thermocline separation. 
Furthermore, measured current shear showed relative current 
in the upper layers to be in a southerly direction on the 
eastern side of the front .and in a westerly direction on 
the western side; currents inside the . boundary region ex- 
hibited diverse directions. 
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Fronts are usually (but not exclusively) determined 
by dynamic factors and are thus convergences (and/or di- 
vergences) of surface currents. They can be caused by 
thermohaline gradients, surface wind systems, or by coastal 
topography, such as coastal fronts occurring, for example, 
along the continental slopes. 

The most common property for recognition of surface 
fronts is sea surface temperature gradient. Shpaykher and 
Moretskiy (1964) have characterized some properties of the 
oceanic Polar front in the following table: 


TABLE 1-1: BASIC PARAMETERS OF OCEANIC FRONTS 


Area - 

Greatest 

Smallest 

Mean 

Width of frontal 
zone, km 

Maximum 

horizontal 

gradient 

Depth of frontal 
divide, m 

Slope angle of 
interface 

wmm 

. 

s, 

% 0 /km 

Arctic 

greatest 

108 

0,39 

0,59 

125 

, ■■ - ■■ 

6' 00" 


smallest 

2 

. 0,01 

0,01 

7 

O' 07" 

56cL 

mean 

21 

. 0,10 . 

0,12 

21 

1 ' 17 " 

Greenland 

greatest 

75 

0,30 

0,08 

1000 

1 ° 08 '16" 


smallest 

20 

0,01 

0,02 

150 

18 ’ 20" 

jca 

mean 

46 

0,11 

0,04 

531 

32 ’42" 

Norwegian 

greatest 

200 

: 0,16 

- 

1500 

1 0 08 ' 4 5" 

Sea 

smallest 

24 

0,02 

- 

100 

12' 44" 

mean 

80 

j ' 

0,07 


745 

3 1 1 19 " 
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As there are seasonal changes in surface temperature, 
frontal discontinuities disappear in some seasons (see 
e . g . , Colton, et. al. , 1975). The coastal fronts are only 
recognizable by temperature gradients (outside upweiling 
areas).- during summer (when coastal waters are warmer than 
offshore waters) and during winter (when coastal waters 'are 
colder than offshore waters). During spring and autumn , -tem- 
perature gradients are frequently masked at coastal fronts. 

The second indicator for surface fronts is water color, caused 
either by different plankton content (quantity of standing crop 
of phytoplankton) or by the amount and type of minercgen- sus- 
pended matter. Salinity discontinuities as well as current 
discontinuities at fronts can be recognized only by special 
measurements and not by remote sensing. Major oceanic fronts 
can be recognized also by changes of waves (surface roughness) . 

Contrary to the transient nature, of meteorological 
fronts, oceanic fronts are relatively stable in respect to 
their mean positions. They frequently represent the time- 
averaged signature of atmospheric features. However, Some 
seasonal shifts, meanaerings and quasi-synoptic changes in 
frontal positions occur (see further Section 3). 

One of the main objectives of the, present study is to 


construct numerical, models for the study of the movement as . 


well as changes in sharpness (gradients) of the fronts as 
affected by. surface driving forces, such as winds. One of 
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the practical reasons for this is that surface- thermal 
fronts can be recognized in cloudless areas from infrared 
satellite observations. However, during the periods when 
the frontal areas under consideration are cloud covered, it 
might be necessary to complement the satellite observations 
with numerical computations on the possible changes in 
frontal positions and sharpness-. 

A further objective of this -study is to model the 
effects of oceanic fronts on surface roughness (waves, wave 
patterns) so that the microwave radiometer gbservations 
from satellites could be correctly interpreted as to which 
surface roughness (wave pattern) changes are due to true 
sharp discontinuities in surface .wind patterns, (meteoro- 
logical fronts) and which are due to major oceanic fronts 
with large current shears. 
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2.0 APPLICATIONS 

Oceanic thermal fronts often coincide;, with, the bound- 
aries of commercially important fish species. . There may 
also be an aggregation of fish near or at the oceanic 'thermal 
fronts (see e.g., Lee, 1952; Hela and Laevastu, 1962). Thus 
it is of economic importance to fisheries to analyze .and pre- 
dict oceanic fronts. 

As the currents of major oceanic fronts, such as the 
Gulf Stream front,, can be strong and opposite in direction, . 
use can. be made of . the currents in navigation for time and 
fuel saving in shipping. 

The vertical thermal and salinity structure (the: tem- 
perature and salinity) change with depth on different sides 
of the oceanic fronts. Thus, the propagation patterns of 
sound will be different on different sides of oceanic f rones 
and the fronts will greatly affect sonar ranges through them 
(Laevastu, Wolff, and LaFond, 1970).' LaFond and LaFond (1971) 
for example, found that temperature inversions at the Cali- 
fornia Front (which is a "medium intensity" front) produce 
sound channels, above the level of increased temperature. Some 
further effects of oceanic fronts' (Polar Front in the Denmark 
Strait) on the propagation of sound have been. described by 
Nunn (1968). Consequently, the knowledge of the positions 


2-1 



of oceanic fronts is of utility to the Navy, for submarine 
as well as anti-submarine warfare. 

Finally, a large part of the energy for the- atmosphere 
originates from the oceans (sensible and latent heat ex- 
change) . This feedback is largest where the difference 
between the properties of the sea surface (temperature .and 
water vapor pressure) and the surface air are large. large 
differences of these properties occur where the tempefauure 
gradients of sea surface temperature are sharp (i.e., at 
oceanic fronts) , (Laevastu and Hamilton, 1972) . Cyclogenesis 
often occurs in the vicinity of major oceanic fronts , and the 
climatological mean position of Arctic and Polar atmospheric 
fronts over the oceans coincides closely with' the correspond- 
ing major oceanic fronts. 
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3.0 SYNOPTIC AND SEASONAL BEHAVIOR 


Although general knowledge indicates that oceanic fronts 
are relatively stationary over synoptic time periods (up to 
3 days) , no detailed synoptic study on the dynamics of off- 
shore oceanic fronts has been carried out. Available studies 
of offshore oceanic fronts using multiple ship traverses show 
that the major fronts are relatively stationary, but include 
relatively large-scale meanders, whereas the minor fronts 
(e.g.. Sargasso Sea fronts) show relatively irregular -and 
small scale meandering and considerable changes of position 
and intensity. 

Colton, et al, 1975, found that surface temperature 
fronts in the Sargasso Sea divide it into a cooler, more 
productive northern part and warmer, less productive southern 
part. They also found a sharp faunal change across this 
front. This surface thermal front in the: Sargasso Sea, 
however, cannot be found during the warm season from July 
to September. One of the main objectives of this study is 
to model and investigate the possible changes of offshore 
oceanic fronts as affected by surface driving forces (mainly 
surface winds) . 

Recent studies of coastal fronts (oceanic fronts near 
continental slopes) show that these fronts are greatly bound 
to bathymetric features (continental slope) , but that changes 
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(small meanderings and mixing) at these fronts occur in a 
tidal rhythm as well as through the action of surface wind.' 

Seasonal migrations of major surface fronts have been 
studied by Laevastu and LaFond (1970) and others. Stefanson 
(1969) found that the near-shore end of a frontal zone 'south- 
east of Iceland appears to move north-south with the season, 
and that it might be displaced as much as '50 'nautical ‘miles . 
Some minor fronts in low latitudes can migrate in some seasons 
as much as 230 miles with an average speed of 7 miles per day 
(He la and Laevastu, 1962). Rossov and Kislyakov (1969) found 
that seasonal and year-to-year changes of the position of the 
Polar Front in the Atlantic can reach 200 nautical miles.. 

Voronina (19 6 2) found that the Antarctic fronts is 
largely determined by surface winds and that it can change 
from a convergence to a divergence if surface winds change 
in the same sense.. Besides changes in posi tion , some sea- 
sonal changes of sharpness (gradients) of the fronts also 
occur . It can be assumed that changes of current in the sur- 
face layers plus changes in surface wind regimes and heat 
exchange processes are responsible for these seasonal changes. 
Considerable seasonal change can occur in the color and 
turbidity of different water types due to seasonal cycles in 
phytoplankton production . 

It is known from empirical observations as well as 
from theoretical considerations that there is considerable 
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variation in the therraocline in frontal zones. These changes 
are largely determined fay the dynamics of the upper, layers 
of the ocean. No detailed studies, nor numerical modeling, 
of the thermpcline changes at frontal regions are, available 
at present. The models (HN models) programmed within ->this 
project will contribute to our knowledge in this area . 

It is well known to seamen that currents 'Strongly affect 
wave characteristics. : If the current runs against, the 'Sea, 
waves are steeper, the length appears shorter and the seas 
are very choppy. Again > if the current runs with the sea, 
the seas are smoother and the waves rounder and longer. Those 
effects can easily be observed in narrows where strong tidal 
currents occur and in channels between islands. Quantitative 
observations at sea are 'scarce on the effects of currents on 
wave shape, length and height. Groen and Dorrestein (1958) 
have constructed graphs which give the ratios of wave length, 
speed, height and steepness of waves in currents with different 
speeds to the same elements in non-flowing water. 

There is a mass transport (wave current) connected with 
waves. Therefore, a long swell (which is causing a "wave 
current") affects also the shorter superimposed waves (sea) 
the same way as currents do. Any swell trains which run 
nearly in opposite direction to the sea can be expected to 
cause steeper, shorter waves and , a swell running in the same 
direction as the sea should cause the sea (wind waves) to be 
less steep, round and longer. This effect has been. treated 
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theoretically by Longue t-Higgins and Stewart (1960, 1961). 
These effects could be expected to occur at oceanic frontal 
regions. Another set of models (spectro-angular wave model 
and wave refraction model) programmed and/or adapted within 
this project will help to evaluate these .effects -quantita- 
tively. 



4.0 NUMERICAL STUDY OF FRONTAL DYNAMICS 

This study is a preliminary; demonstration of the ap- 
plicability of a hydrodynamical-numerical model to the simu- 
lation of oceanic f rontogensis . As the model has been 
described elsewhere in this report, this .section is devoted 
to a description of three initial test (siTanintion)--riins>- and 
a final test run for a case with unusually strong initial 
thermal gradient and convergent winds. 

In all cases a two layer ocean was simulated in which 
the lower layer was initially at rest. In each case a current 
field, initially specified in the top layer, developed in 
both the top and bottom layers in response to an imposed 
wind field and to the initial imbalance in the mass and 
velocity fields. As the initial state of the pressure (depth) 
and density structure needed to balance the specified motion 
fields are not known, the model must be run for a sufficient 
number of real-time hours to establish proper dynamic balance 
between mass and motion fields.- In all four test cases 
described here, run times were extended to the point where 
balance .was achieved. This problem and the way in which 
balance is finally realized within the model itself is common 
to all Primitive Equation models (including those used for 
routine atmospheric prediction) . 
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In order to demonstrate the use of the HN model for rou- 
tine conditions of initial current structure and surface wind 
stress, three different simulation runs were conducted in an 
initial test series. These consisted of a twelve-hour simu- 
lation of a current shear front (Case I) , an eighteen-hour 
simulation of a surface wind convergence front.. (Case II) , ...and 
a twleve-hour simulation of a cyclonic wind field moving 
across a current shear front (Case III).. The several impor- 
tant physical constraints which differed in the simulation 
runs are listed in Table 4-1 . 

In Case I the initial current and steady state wind 
field, uniform in the east-west direction, are shown in 
Figure 4-1 . The initial current was zonal, non-diver gent 
and had linear shear in the northern and southern portions 
of the grid. The wind field, also specified as a. shear flow, 
was divergent south of the front and convergent north of 
it. It should be noted that at the front the initial cur- 
rent and wind field were set to zero. 

In Case II the initial current field was prescribed 
as uniformly flowing to the west at 2.5 cm sec ^ . The wind 
field, shown in Figure 4-2, simulated a convergent wind 
stress front. North of the front, the winds were divergent.; 
and south of the front, they were convergent. In Case II, 
as in Case. I, the wind field varied only in the north-south 
direction . 
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TABLE 4-1: 


PARAMETER SPECIFICATIONS 


PARAMETER 

CASE I 

CASE II; III 

TOTAL DEPTH 

300 m. 

300 m. 

GRID MESH 

10 km. 

10 km. 

TIME STEP 

60 sec. 

60 sec. 

THICKNESS OF SURFACE 
LAYER ..... 

50 m. 

75 m. 

DENSITY OF SURFACE , 
LAYER 

1.026 g cm 2 

1.025 g cm 2 

DENSITY OF BOTTOM 

„ „ -3 

-3 

LAYER 

1.029 g cm 

1.027 c cm 

LATITUDE OF CORIOLIS 
FORCE 

30° N 

45° N 

SURFACE LAYER SMOOTHER ( 15 

.998 

.998 

BOTTOM LAYER SMOOTHER 

.920 

.920 

EDDY DIFFUSION 
COEFFICIENT (2) 

4 

5x10 

4 

5x10 


(1) These parameters simulate eddy viscosity by considerin 
the velocities at surrounding grid points. 

(2) The value of this parameter in Case III was 5x10^. 
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WIND FIELD 
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FIGURE 4- 


1 m 


-1 

sec 


CURRENT FIELD 

1 cm sec 


. : WIND AND CURRENT . FIENDS FOR CASE . I 





In Case III the initial current field was that used 
in Case I. An initial uniform wind field of 5.4 m sec - ^ from 
112° (south of east) was prescribed over the grid. As the 
simulation proceeded this wind field was modified by a cyclon- 
ic disturbance. This disturbance was of a radius of 60 nau- 
tical miles and moved from the southern boundary to the North- 
west. 

The discussion of the results of these inital cases is 
limited to a description of the dynamic behavior of the two 
layer ocean-. In Case I a shear- front , in both the wind field 
and current field was prescribed. The front, as defined by 
the slope of the interface between the top and bottom layers, 
is shown in Figure 4-3. The top layer deepens under the influ- 
ence of the convergent wind north of the front and shoals 
under the divergent regime south of the front. The initial 
trough in the interface occurring at the location of zero 
velocity in the Wind and current fields becomes less pronounced 
by the end of the twelve hour simulation run. 

The currents in the surface layer in Case I at locations 
north and south of the front are shown in Figure 4-4 in- the 
form of progressive vector diagrams. From this figure,, if 
can be seen that the primary mode of the: response of the sim- 
ulated ocean to the initial conditions and the steady-state 
wind field is that of a rotary current. 
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Deviation of Interface (CM) 
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N = 1 
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15 
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FIGURE 4-3: 


DEVIATION OF INTERFACE FROM INITIAL POSITION 
ALONG M = 15 CASE I 
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12 Hr 



CURRENTS NORTH 
OF FRONT 


1 cm sec 

- -1 
5 cm sec 




CURRENTS 


SOUTH OF 


FRONT 


FIGURE 4-4: CURRENTS IN TOP LAYER - CASE I. 



Case II simulated a convergence front, the development 
of which is shown in Figure 4-5. In this run the front 
develops without a trough in the interface as happened in 
Case I. This fact can. be attributed to the lack of a null 
current line in the initial current field in Case II. 

After twelve hours the slope of the interface in this run 
was about 50% of that in Case I. Also, it appears that 
the top of the front is eroding northward from its initial 
position resulting in a gradually, decreasing frontal slope. 

The currents north and south of the front in Case IX 
are again depicted by progressive vector diagrams and are 
shown in Figure 4-6. While the rotary current response is 
the same as seen in Case I, there is not the initial change 
in direction across the front in Case II due to the uniform 
initial currents' in this case. Again no change was found 
in the initial temperature distribution in Case II . 

Case III simulated the interaction of a time dependent . 
windfield and an initial shear current. The development of 
the front is seen in Figure 4-7. In this case the front 
develops through the first four hours as in Case I- At 
six hours however, the topography of the interface becomes 
distinctly different from that of Case I. This interface 
topography .continues to evolve differently through, the. remain- 
ing six hours of the run. The general form of the front is 
similar to, that of Case I . in which the interface deepens 
in the north and shoals to the south. 
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Deviation of Interface (CM) 



FIGURE 4-5: DEVIATION OF INTERFACE FROM INITIAL 

POSITION ALONG M = 15 CASE II 
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FIGURE 4-6 : CURRENTS IN TOP LAYER 

CASE II 
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Deviation of Interface (CM) 



FIGURE 4-7: DEVIATION OF INTERFACE FROM. 

INITIAL POSITION ALONG M=15 
CASE III. 

4-12 ' ' ; 


The currents in the surface layer in the northern and 
southern regions are shown in Figure 4-8. Because in 
Case III the wind field is not uniform in the east-west 
direction, the currents are not as representative of the 
region concerned as they are for Cases I and II. By com- 
paring Figures 4-4 and 4-8, the impact of the different 
wind regimes can be seen. In the northern region the rotary 
current is rather strongly deformed compared to Case I due 
to the differing wind regimes. 

Since the initial test runs did not show appreciable 
changes in temperature structure in the upper layer (because 
of the weak initial thermal gradient and . ;comp aba ti. ve 1 y ':;sbott 
run time), an additional test case was formulated to demon- 
strate that temperature changes would occur. In this case, 
a strong front was specified initially and current/wind 
forces were specified which should cause a gradual weakening 
of the frontal intensity --especially in the frontal fringes. 

Figure 4-9 shows the initial north-south temperature 
profile (homogeneous cold and warm water masses with a sharp 
gradient in between) and the change in surface layer temperature 
Which occurred after 120 hours of run time. Notice that some 
blurring can be seen at the cold edge but the largest changes 
(almost 0.6°C) occurred at the warm edge of the front. This 
non-symetrical pattern is due to the way. the initial winds 
and currents were defined. Figures 4-10 and 4-11 show the 
initial and final temperature values. 
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FIGURE 4-8: CURRENT IN TOP LAYER 

. CASE III 
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FIGURE 4-9: ORIGINAL TEMPERATURE PROFILE AND 

TEMPERATURE C FLANGE AFTER 120 HRS. 
CASE IV. 
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FIGURE 4-10. INITIAL TEMPERATURE FIELD 


4-16 




5 . 0 NUMERICAL STUDIES OF SEA STATE CHANGES 

One basic objective of this overall,;.studY:wa.s to ^deter- 
mine the influence of oceanic fronts upon sea state - es- 
pecially the effect of fronts on wave energy distributions 
and/or variations in wave steepness which might be detected 
by satellites. Two approaches to this problem were utilized: 
(a) a spectral wave model (see App . 3 ) was. used to study 
sea-state energy distributions around - an atmospheric con- 
vergence zone which would support/amplify .an oceanic front 
and (b) a model described by Johnson (1947) was used to study 
wave changes caused by intersection with a current boundary, 
i . e . , a numerical wave refraction model was used as. described 
in App . C . 

5 . 1 Wind Wave Changes Across Oceanic Fronts 

a. Experiment Design 

The French D.S.A. V spec troangular wave; model was 
used to generate typical wind wave fields resulting from an 
atmospheric convergence situation and, also, the passage of 
a tropical cyclone. The description of this numerical model 
is found in App. B of this report,. To reproduce the wind 
fields, the 29 row x 33 column grid was used, a wind speed 
of 20 kts from the ESE occurring at row 14 and above and a 
wind speed of 10 kts from the ESE at row 15 and below. For 
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the tropical cyclone experiment, a cyclone of radius 200 nmi 
and maximum wind speeds of 60 kts , at the wall of the eye 
was prescribed. The wind speed profile was based on: 




= V 


max 


V - V 
max Back 

In ( R /R . ) 

eye circ 


In (D/R ) 
eye 


where: 


V 

Back 

R 

eye 

R . 
circ 

V 

max 


is the velocity at distance D 

is the velocity of the background flow 

is the radius at the eye 

is the ■radius of the storm 

is the maximum wind velocity 


In both experiments the grid mesh. length was set at 
40 nmi. The time step between introduction of wind fields 
was 3 hr. The program was run out through 40 hours of com- 
putations in both instances . 

b. Case 1 : Convergent. Wind Field Results 

Figures 5-1, 5-2 and 5-3 depict a subsegment across 
the front of the fields for direction of maximum 

energy, and period of maximum energy, after 33 hours of 
computation. Higher wave heights are evident on the northern 
side of. the front as well as longer periods. The fact that 
the heights decrease, 10 ft. to 7 ft. , and the periods shorten, 
8 sec. to 6 see., as the convergence is approached is due to 
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FIGURE 5-1 

H 1/ , 1q Wave Height (ft) at 
Hour 33 of Convergent Case 
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FIGURE 5-2 

Wave Direction of Maximum Energy at 
Hour 33 of Convergent Case 
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FIGURE 5-3 

Period of Maximum Energy (sec) at 
Hour 33 of Convergent Case • . 
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the combined effects of increased wave-interference and 
dampening. On the southern side, the lighter wind speeds 
result in a comparatively flat H 2./] 0 °f ^ ft. w ith 

a period of maximum energy of 4 sec. 

The direction field shows two bands of anomalous di- 
rections of maximum energy. These are clarified when the 
spectral information at the three special points are in- 
spected on Figure 5-4. As can be seen , . the direction of 
maximum energy , when .. moving from points 1 to 3 , undergoes 
a dampening of ENE components, while- the NE components remain 
large enough to become the direction of maximum energy/ The 
overall pattern shows a decrease in energy when crossing the 
front accompanied by a corresponding shift in the direction. 
Likewise, the period of maximum' energy shortens once the 
energy decreases while propagating southward across the 
front . 

c . Case 2 : Tropical -Cyclone .Results . 

For this, experiment results for three different 
times are presented (21, 27 and 33 hours of computation). 
These times correspond to the passage of the storm -across 
the convergent front described in the /previous section.. . 

In this case, the storm was initiated at t=0 and allowed 
to travel across the wave field, produced by the convergent 
wind fields. 
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FIGURE 5-4 : SPECTRAL COMPONENTS OF SP 
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Figures 5-5 through 5-7 depict the surface wave para- 
meter fields generated by the wind field shown, in Figure 
5-8. As can be seen, much higher values of occur 

around the storm and. to its right where the storm velocity 
and background flow coincide. The directions of maximum 
energy take on a cyclonic appearance. (Values of : "382° " 
mean that no single direction has a predominance of energy.) 
The same pattern is evident in the wave period field, with 
higher periods occuring to the right and ahead of the storm. 
Here, points of no predominant period of maximum energy 
are left blank. 

Much the same results are evident from the results at 
27 hours, Figures 5-9 through 5-12. At this time the storm 
is in the process of crossing the convergence zone. From 
Figure 5-9, it can be seen that the values of ^^/1Q ^ ave 
spread across the zone. The directional field still has a 
cyclonic nature and has almost totally obscured the convergent 
initial pattern. In the wave period field the high periods 
are evident moving ahead of the storm path. An interesting 
aspect is the spread of 11 sec. periods along both sides of 
the convergent zone. 

Figures 5-13 through 5-16 depict the results after 33 
hours of computation. : At this time the storm is largely 
across the convergent zone and moving away. As can be 
seen, the wave height field has decreased substantially, 
reverting to the original state. The directional field. 
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FIGURE 5-5 

H 1/i0 Wave Height Field (ft) at 
Hour 21 of Tropical Storm Case 
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FIGURE 5-6 

Direction of Maximum Energy (° from) at 
Hour 21 of Tropical Storm Case 
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FIGURE 5-7 

Period of Maximum Energy (sec) at 
Hour 21 of Tropical Storm Case 
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FIGURE 5-8 Wind Speed (m sec ) and Direction ( 0 from) for 

Hour 21 of Tropical Storm Case 
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FIGURE 5-9 

H l/10 Wave Field (H) at 

Hour 27 of Tropical Storm Case 
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FIGURE 5-10 

Direction of Maximum Energy ( ° from) - at 
Hour 27 of Tropical Storm Case 
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FIGURE 5-11 

Period of Maximum Energy (sec) at 
Hour 2 7 of Tropical Storm Case 
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IGURE 5-12 Wind Speed (m sec ) and Direction t° from) for 

Hour 2:7 of the Tropical Storm Case 
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FIGURE 5-14 

Direction of Maximum Energy (° from) at 
Hour 33 of Tropical Storm Case 
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FIGURE 5-15 

Period of Maximum Energy (sec) at 
Hour 33 of Tropical Storm Else 
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FIGURE 5-16 Wind Speed (m sec ) and Direction ( from) for 

Hour 33 of the Tropical Storm Case 



however, is still quite extensively altered due to the passage 
of the storm. The periods of maximum energy have begun to 
revert to their original state with the. longer periods occur - 
ing at the convergence zone. 

Figures 5-17 thorough 5-19 depict the directional and 
periodic energy spectra for the three special points as 
functions of time. Despite their locations which account for 
initial energy differences, all three show the same overall 
patterns in the energy versus direction spectra. A shift 
in direction to the southeast is evident in all three and, 
except for the northern^-most point, all experience a decrease 
in energy. In the case of period of maximum energy. Point 1 
experiences an increase to longer periods, Points 2 and .3 a 
decrease. Only Point 1 experiences an increase in energy in 
either spectra. 

5.2 The Effects of Oceanic Fronts on Wave Refraction . 

The theory presented in this section is based upon 
Johnson's (1947) methods for calculating the effects at 
oceanic velocity fronts on the refraction of waves. The 
basic equations are presented in a following section (App. C) 
in which the full computer program is described. This section 
will deal only with the results of Johnson's theoretical 
study and the numerical results based upon his work. 
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The core of Johnson's results are summarized by 
Figures 5-20 through 5-22. Given m, the ratio of current 
velocity over the wave ray speed, Figure 5-20 shows that the 
incident angle, a, is transformed into angle B dependent 
upon this ratio. The largest alterations occuring in cases 
where the two are of equal magnitude, with the current 
either opposing (negative on) or conciding (positive on) 
with the ray direction. 

Figure 5-21 shows the effect of the currents upon the 
steepness of the waves. Here the largest effects are found 
at m=-l. This value is relatively .uncommon in deep water 
cases but can easily be encountered in near-coastal areas 
such as river outlets. As can be seen from Figure 5-22, the 
maximum effect upon the wavelength occurred at .im=+l, m=-l 
resulting in a shortening of the wavelength. From the above 
results, it is clear that wave height is increased with ne- 
gative values of m and median values of 

Two examples of the results of the numerical model 
are shown in Figures 5-23 and 5-24. Table 5-1 summarizes 
the results of all computer -runs . From this it can be 
seen that the maximum height change occurred at a=30, 
m=-.064. The maximum directional shift was —9° at a=60, 
m=-.064. The maximum steepness occurred at a=45, m=- . 064 . 

It can then be concluded that the most severe alterations 
at wave fronts occur when large opposing currents are 
encountered. 
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FIGURE 5-20 : THE EFFECT. OF WAVE DIRECTION ON 

THE REFRACTION OF WAVES BY 
CURRENTS 
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FIGURE 5-21: . STEEPNESS FACTOR IN REFRACTION 
OF WAVES BY CURRENTS •" 
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FIGURE 5-22: WAVE LENGTH CHANGES DUE TO 

REFRACTION OF WAVES BY 
CURRENTS 
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FIGURE 5-23 : EXAMPLE OF COMPUTER OUTPUT, a=45° , C=-50 cm/sec 
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FIGURE 5-24 s EXAMPLE OF COMPUTER OUTPUT, a» 45 °> C «-50 cnt/aoc 





6.° THE EFFECTS OF' OCEANIC FRONTS ON, AIR/SEA EXCHANGE 

PROCESSES 

It is known that the ocean-atmosphere heat and moisture 
exchange (feedback) is different on opposite sides of oceanic 
(and atmospheric) fronts (e . g . , Laevastu and Hamilton, 1972). 
This relatively abrupt spatial exchange of heat and mois- 
ture (evaporation) might affect various .observations and 
sensings by satellites of the near-surface atmosphere as well 
as the sea surface. 

In order to study and demonstrate quantitatively the 
feedback of energy from the oceans to the atmosphere, the 
best exchange formulas, verified in . various .past research, 
were programmed into an integrated computer program (App. D 
of this report) . In this program, sea surface temperature 
(T w ) and winds (V) were prescribed; all other parameters 
were computed. The results of a 2 4 -hour "steady state" com- 
putation are shown in Figure 6-1. 

The heat and moisture feedback from the ocean depends 
greatly on the removal of these quantities through turbu- 
lent transfer, -which is a function of surface wind speed. 
Furthermore , the quantitative turbulent transfer is also 
dependent on the difference between the sea surface tempera- 
ture and water vapor pressure and the corresponding proper- 
ties of the near-surface air, (T T and e - e- ) . 

w a w a 

6r-l 



The properties of the surface air (temperature and water 
vapor pressure) adjust relatively rapidly to the corresponding 
properties of the sea surface, so that an "equilibrium dif- 
ference" is established within five to six hours. The adjust- 
ment of the properties of surface air to the corresponding 
properties of the sea surface, and the resulting "equilibrium 

difference" is dependent on the rate of change of sea surface 

3T 

w 

temperature under the trajectory of surface air (V t — -).. 

r dr 

Thus, the surface wind speed and its direction in relation 
to the sea. surface isotherms determines the heat and mois- 
ture exchange (see Figure 6-1) ; the sea surface temperature 
gradient influences this exchange to a lesser degree. 

Consequently, if the change of surface winds (both 
with respect to speed and direction) closely coincides 
with the oceanic front, a relatively sharp gradient of latent 
(Q q ) and sensible (Q^) heat exchange occurs at the front 
(Figure 6-1). However, a comparable gradient of Q g and 
can also be found at atmospheric fronts which do not coincide 
with oceanic fronts. 
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7.0 SUMMARY 


A two-layer Hydrodynamic Numerical (HN) model of the 
Hansen (primitive equation) type was used to show how ocean 
temperature structure . and near-surface currents are .'influ- 
enced by surface wind stresses in a case where a convergent 
wind field coincided with an ocean front. The ;;Erenr*.h . -S.pect.ro 
angular Wave Model was used to show differences in spectral 
wave energy across an ocean front being (a) maintained by 
surface wind convergence, and (b) disturbed by a moving 
cyclone. The studies of Johnson (1947) on wave refraction 
due to currents were programmed to demonstrate chah'ge'S >wh-ich 
can be expected in wave direction, height, and steepness 
due to intersection with a strong current, stream (front). 
Finally, a computer program was developed to show changes 
in transport of latent and sensible heat due to air/sea 
temperature. differences across an ocean front. 

These studies indicate that surface properties vary 
sufficiently in the vicinity of major fronts to permit^de- 
tection by satellite. 
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APPENDIX A 


HYDRODYNAMICAL-NUMERICAL (HN) MODEL 
FOR THE- STUDY OF OCEANIC FRONTS 

CONTENTS 

1. Basic Hydrodynamical Equations' arid Their 
Finite Difference Forms 

2. Model Inputs and Other Parameters 
and Program Organization 

3. Symbols Used in the Program 
Program listing separately submitted 



1 . 


Basic Hydrodynamical Equations and Their Finite 
Difference Forms * 


A two-layer Hydrodynamical -Numerical (HN) model of 
Walter Hansen's type was used in the present study. The 
individual layers are vertically integrated. The basic 
equations are as follows: 


+ H , (Ul ) + H ' (VI ) = 0 
^1 ^2 ul x vl y 

C- + H _ (U2 ) + H - (V2 ) = 0 
^2 u2 x v2 y 


Ul + 


r Tui 2 + VI 2 


H 


Ul 


Ul - f VI — g c lx = K (x) 


U2 + 


r >/U2 2 + V2 2 


H 


u2 


U2 - fV2 - g £| ?lx - g(l-£|)C 2x 


VI + 


r /ul 2 +V1 2 


H 


vl 


VI + fUl - g 5 = K (y) 


V2 + 


r t/U2 2 + V2 2 


H 


v2 


pl 


... pl^ 


v2 + fu 2 - g ^ c ly - ? 2y 


where : 

- surface elevation; 

£ 2 - deviation of MLD - (mixed layer depth) from mean 
U1,V1 - u,v components in upper layer 
U2 ,V2 - u,v components in lower layer 
r - friction coefficient 
f - Coriolis parameter 
g - acceleration of gravity 
H - depth 

P 1 ,P 2 - ( ^ ens l t l es 

K tx) , - external forces 
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There are two inter-dependent continuity equations, 
one for each layer. These compute the change of sea level, 
the change of depth and thickness of the layers. 

The equations of motion for each layer are vertically 
integrated through this given layer. The lower layer is 
driven by internal friction and by pressure gradients.. 

The finite difference forms used for computation are: 


-t+T _ -t-T _ Tr„t 

C / _ _ 1 % n t H. 


u 


- PI 


(n ,m, 1) ^ (n,m, 1) Z u (n,m, 1) (n,m,l) u (n.,m-l ,1) 


u 1 ; H fc . T V - H fc , . . 

(n,m-l , 1) v (n-1 ,m, 1 ) (n-l,m,l) v(n,m,,l) 


v 1 ; ~ { H 1 , - H t / . 

(n,m,l) Z u(n,m,2) (n,m,2) u.(n ,m-l , 2 ) 


t t t t 

U (n,m-1,2) + H v (n-1 ,m, 2) V (n-l,m,.2) " H v.(n,m,2) 


V (n,m,2)> 


t+T ■ ' -t-T _ T_ r t . TJ t. _ H t 

^ (n,m, 1) ^ (n ,m, 2 } Z u ( n , m , 2 ) (n,m,2) u (n,m-l,2) 

t t t t 

U (n,m-l , 2 ) + H v(n-l,m,2) + V (n-l,m,2) " H v(n,m,2) 


V, } 

(n,m, 2) 




T t+2T 


T t+2T - 


- 1 *t 

U (n,m,l) + 2 fxV (n,m,l) 

? 1 7 +X \\ } + 2tX 1 7 +2t . 

s (n,m,l) (n,m) 


I£ fr t+T 
Z ' 4 (n ,m+l , 1 ) 
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^n, m ,2, 2 + V (n,m,2) 2} 


,t + 2T 
(n,m, 2) 


C (n,m,2) + 2xf V (n,m,2) £ { p2 Cc (n,m+l.,l) 


_ H f£i r r t+T 


r t+x n _ 13. r ri - fill rr t+T 

? (n,m,l) Jx £ 1 L p2 J u (n,m+l,2) 


_ C t+T 1} 

5 (n,m,2) u 


V (n,m,l) ^ t 2xr / H v (n,m,l) ^ V (n,m,l) +U (n,m,l) } 


o-rf r :** 1 - IS. f r t+T r t+T > 

2xf U (n,m / 1) £ 1 Mn,m,l) ^ (n+1 ,m , 1 ) i 


+2xY 


t+2x 

(n,m) 


t+2t _ ri 1 J ^ _l_ TJ * 

V (n ,m, 2) {1 “ [2x /H v(n,m,2) ] V (n.,m,2) U (n,m,2) 


2 1 


o -F n* t -IS. r£i rr t+x _ r t+T t 1 

~ 2xf U (n,m, 2) £ { p2 U (n,m,l) ; (n+l,m, 1). J 


^ [? t+T 


t+T 


z L L ~ 'p2 J lt3 (n,m,2) C (n+1 ,m, 2 ) - 1 } 

Th.e computation of U and U is done as In the single-la.ye 


model 


U ^n,m) ~ a °^n,m) + 4 {U (n-l,m) + U (n+l,m) + U (n,m+1) 


+ U, , v ) 

(n ,m-l) 


U (n,m) 4 ^ U (n,m-1) + U (n+l,m-l) + U (n,m) (n+l,m) 


A-3 



The computation of V and V are analogous to the correspond 
ing U computations above. The actual depth (and thickness of 


the layer) is computed with the following equations: 


H t+2x , h , 1 r t+x _ t+ t , 

u (n ,m, 1) u(n,m,l) 2 1 ^(n,m,l) ^(n,m+l,l)-‘ 


H t+2x _ h .If _t+x t+x x 

v(n,m,l) n v(n,m,l) 2 1 ? (n,m,l) ^ (n+1 ,m, 1) ' 

The following symbols were used in the finite di-fe-Ierence 


formulas above 

x »y 

t 

U,V 
u , v 


h 


H 



X,Y 

9 

f 


space coordinates 
time 

components of velocity 

indicators of U and V point (location) in 
the grid 

initial depth (when £ = 0) 

surface elevation (for the second layer, it 
indicates the deviation of the depth of the 
layer from its prescribed mean value) 

total depth (H = h + 5) 

depths at u and v points respectively 

components.^ of external forces 
acceleration of gravity 
Coriolis parameter 


r friction coefficient (bottom stress for lower 

layer) 

a coefficient of . horizontal eddy viscosity 

(also acts as a smoothing coefficient) 

n ,m coordinates of the grid point, 1 and 2 

indicate the first (surface) or second 
(deep) layer 

x half time step 


% half grid length 

Pf,p2 densities of first and second layer 
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As the HN model provides a current vector at each grid 
point at each time step, one can solve the interacting -ad- 
vection and diffusion problem using a time - stepping : f i n-i be 


difference method. In the frontal model the distribution of 
surface layer temperature will be advected and diffused to 
obtain the changes in position and intensity of the ■surface 
front as affected by wind-induced current components,, 'ijhe 


following restraints and considerations apply to the advection/ 


diffusion computations: 

(a) Vertical diffusion between layers is neglected 

(b) A Fickian type diffusion with basic Lagrangian 
approach, and with constant diffusivity is used 
since the advection and diffusion equations are 
solved separately in short, interlocking time-steps. 


(c) The Austausch (eddy diffusion) coefficient is 
related to time-step and grid size used. The 
advection is computed linearly in the following 
finite difference form: 


s t 

- X 

u t+T 

(s t ; ; 

- S 1 ) 

- x 

v t+x 

n,m 


n,m 

\ n,m 

n , m+ 1 J 


n ,m 



i 


where: Sn,m-1 or n,m+l (respectively n-1, n+1) are used, 

depending on the direction (sign) of U and V. 
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The following diffusion formula for temperature h a .s 
been used; heat is conserved provided no decay or heat 
exchange formula is applied and the correct Austausch co- 
efficient is chosen. 


S t+T = S t + ^ (S t . + S t ' 

n,m n,m ^2 n,m ^2 ^ n-l,m n,m-l 


+ S 1 . + S t . - 4S t ) + ~ [S t . 

n+l,m n,m+l n,m J ^ 2 ^n-l,m-l 


+ S n-l,m+l + S (n+l,m-l) + S (n+l,m+l)l 


This scheme requires that 


UAt 

A1 


< 1 


The Austausch coefficient A is taken from the graph 
on Figure:!; 


1000 A 
At (graph) 
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Figure 1 also shows some other horizontal dimension 
(grid size) - Austausch coefficient relations found experi- 
mentally by different researchers. 

The maximum length of the time step in the numerical 
model is determined by the grid length and maximum depth 
according to the Courant-Friedrichs-Lewy criterion: 


0 . 5A t < 


0 . 5 £ 

/VC 


where At is time step (sec); Z is grid length (cm) ; g is 
acceleration of gravity; and H max is the maximum depth in 
the area of computation (cm) . Any attempt to increase the 
time step above this criterion results in "blow up" or 
computational instability. 

The selection of the size of the area and the distance 
between grid points is usually based on requirements concern- 
ing topography, accuracy, and availability of computer core 
memory. A^ 10 km grid size was used in the present model. 

The grid net (staggered grid) is shown in Figure 2. 

It consists of three different sets of grid points: (1) the 

water elevation points (Z) at the intersection of the grid; 

(2) the point for U-velocity component to the right and (3) 
the V-velocity component below the corresponding Z point. . 
Each of these three points hah the same coordinate designa- 
tion (n ,m) . 

As this model has been specially adapted 
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for the frontal study, no rotation of the grid from a normal 
position is permitted. 



FIGURE 2: SCHEME OF THE GRID NET 


The solution to the hydrodynamical equations is 
uniquely determined if all boundaries are either closed or 
prescribed. In the present study, the boundaries must be, 
however, left open for in- and out-flow. Much effort has 
been spent in the past to find theoretically and practically 
(from computational point of view) satisfactory solutions to 
open boundary treatment. Two relatively simple approaches 
have proven most useful . In one approach, the variables 
(U,V and Z) are computed out to one row or column from the 
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boundary and these next-to-boundary values are transferred 
to the boundary at each time step. The second approach, 
which has verified somewhat better, is to compute the 
Variables up to the open boundaries and then replace the 
one missing parameter (which can be either U, V or Z) with 
its. value at the boundary. This method is commonly - re- 
ferred to as the "leaky boundary" method. 

2 . Model Inputs and Other Parameters 

and Program Organization 

The smoothing coefficient (ALPHA in the program) acts 
as a horizontal viscosity coefficient. Normal value for 
ALPHA in surface layers is between 0.985 and 0.99. For 
the second layer a somewhat lower value (0.92) is used to 
suppress some computationally caused, gravity waves . 

The bottom friction is a non-linear (quadratic) quantity 
with a coefficient r of. 0.003. The value for the internal 
friction coefficient is somewhat uncertain. Numerical exper- 
imentations have shown that 0.0015 gives verifiable results. 

If the model is applied in real conditions, depth 
values, as well as initial layer thickness values must be 
digitized at U and V points. A sea/ land table is read in at 
Z points... The coast must pass through U and V points and 
no flow into or out from . the coast is allowed. This condition 
is not pertinent to .the frontal model, as it has four open 
boundaries. 
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In the present model the permanent current (e.g. 
thermohaline current) , the wind and surface temperature 
are prescribed at each grid point. The wind speed is given 
on an input card (or in a data statement in present model) 
in m sec The wind direction must be given in grid coordi- 

nates and in the direction toward which the wind is^hlpwing.. 
Thus, a west wind is coming from -x toward +x and should 
be given an angle of Oi; a north wind has an angle of 270°; 
east wind, 180°, and south xvdnd , 9 0°.: 

The FORTRAN program is organized into 'eight subroutines . 
The common blocks are identical in all subroutines. The main 
program calls the major subroutines and counts the time. Sub- 
routine J02 is for initialization of the fields; it inputs 
and computes a number of derived parameters used in other 
subroutines. It is called only once. 

Subroutine VECTOR is called by subroutine JO 4 and compute 
the. resultant direction of the current from U and V components 
during output times. Subroutine J04 is for output and is 
called at any desired time interval (Tl) . The latter is 
given in the data statement. This subroutine also writes para 
meter values for special (selected) points on a tape for 
output at the end of computations (using subroutine J08) . 

Subroutine JOS is the main computational subroutine 
and is called each time step. This subroutine also calls 
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wind subroutine SOI and the current and temperature 
initialization subroutine (CURRENT and TEMPS) . Subroutine 
SOI calls the wind subroutine (WINDS.) and computes the 
wind current component. 

Subroutine WINDS contains various (different) frontal 
wind generation sections, different types of initial frontal 
surface current generation sections and the initial surface 
temperature distribution section. 

References 

QCEANOGR. DEPT. EPRF, 1974. A vertically integrated Hydro-. 
dynamical-Numerical Model (W. Hansen : type) . 
ENVPREDRSCHFAC Tech. Note 1-74. 

LAEVASTU, T. , 1974. A multi-layer Hydrodynamical-Numerical 

Model (W. Hansen type) . ENVPREDRSCHFAC Tech. Note. 2-74 
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3. 


System Variables 


3 . 1 Description of Input Parameters 

The following parameters are set internally in subrou- 
tine J02 in DATA statements: ; 


AGFN 
* ALPHA 

*C 

Cl 


DEBUG 


DL 


D+ 


F 


G 

HTL1 


NOT USED 

Velocity field smoothing parameter or pseudo 
eddy vescosity for top layer; ALPHA- .993 

Drag coefficient for computation of wind stress 
C = 3.2 x 10~ 6 

Flag to indicate whether or not wind stress is 
imposed; Cl = 1 for wind stress, Cl = 0 for no 
wind stress 

Flag used to print and initial current, wind and 
temperature fields; DEBUG = 1 prints, these fields 
at beginning of the program, DEBUG— 0 suppresses 
the print 

1/2 distance between nodes on the computational 

5 

grid; DL = 10 cm 

1/2 the tine between iterations of the integration 
scheme; DT = 60 sec 

Coriolis parameter; F = 2usin$ 

u = angular rotation of the earth; <j> = latitude; 

(F = 1.0312 x 10 _4 at Card 29 of J02) 

-2 

Acceleration due to gravity; G = 970. cm sec 

Initial depth, (thickness) of the top layer; 

HTLl = 7500.0 cm 


IUE NOT USED 

IZE NOT USED 

JA Should be zero always 


KKE 


NOT USED 
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LI 


Reset at Card 341 of J05 to LI = 1 as temperature 
advection is calculated (may be removed from data 
statement) 

ME Number of columns in the computational grid; 

ME = 30 

NG NOT USED 

NE Number of rows in the computational grid; 

NE = 25 

NS NOT USED 


NURU 

*R 

RBETA 

ROL 


DUMMY variable; leave at NURU - 6 
Bottom friction coefficient,-; R = 3.x. ,10 


-3 


Used only in PRINT statement in subroutine J02 
(card 63) which may be removed 


NOT USED 


SI Time (sec) at which wind field 'is set to zero; 

4 

SI = 36 x 10 (SI must be greater than TE for 
wind input throughout run) 

SIGMA NOT USED 

T Time index - always set at zero, ■initially 

TE Total time for integration of equations 

(Simulated real time) ; TE = 43200. sec 


TW Time (sec) at which wind field is first imposed; 

TW = 100.0 sec 

Tl Output (printer) interval in seconds; T == 3600 

implies output parameters will be printed at 
1-hour intervals 

T2 Leave at zero always 

T3 NOT USED 
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The following are set in DATA statements in J05: 

*ALP Velocity field smoothing parameter or pseudo 

eddy, - viscosity for bottom layer.; ALP = 0.92 

IFLG Flag used to initialize current filed on first 

iteration; IFLG ^ 0 results in bypassing this 
initialization 

IFLG2 Flag used to set up temperature field on first 

iteration; IFLG2 = 0 results in bypassing the 
temperature field initialization 

*RHOl Density of the top layer; RHO = 1.025 g cm ^ 

*RH02 Density of the lower layer; RH02 = 1.027 g cm ^ 

The following- is defined in J05 at card 350 : 

AUS Eddy diffusion coefficient; AUS = 5.0 x 10 

*The range of these variables is given in Table A-l. 
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3 . 2 Initialization Options 


The following are used to select an option in the ini- 
tialization of the current, temperature and wind fields. 


Card 109 selects current (cm sec field 
initialization 

Currents parallel to the front, but in opposing 
directions from one side of the 'front to the 
other. (Case I and III Current Field) 

Currents converging and decreasing toward front 
(not used to date) 

Uniform east west currents at 2.5 cm sec ^ (Case II ) 

Currents converging and increasing toward front 
(not used to date) 

Currents diverging from front (not used to date) 


In J05 : 

Card 332 is used to select the temperature field 
initialization 

KEY=1 This is the only option at present for the initiali- 

zation of the temperature field. The temperature 
field is homogenous in the east- west direction . 1 


In J01 : 

Card 20 is used to select the wind field initialization 

KEY-1 Case I wind field divergent south .of front, 

convergent north of front 

KEY=2 (has not been used to date) 

KEY=3 Call to subroutine HURRW for "the generation of a 

moving wind field 


In J05 : 

KEY=1 

KEY=2 

KEY=3 

KEY=4 

KEY=5 
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3.3 Description of Output Parameters 

The normal output (DEBUG=0) consists of initial 
printout of the values of : DT, DL, ALPHA, R, : F , 

SIGMA, RBETA, LAMDA. 

SIGMA (NOT USED) 

RBETA (NOT USED) 

LAMDA is the variable C described in 3 . l,,of this 

appendix 

Following the printout of the above, these fields aire 
printed at the interval specified by T1 : 

Z . Height of the free surface above the initial 

position in CM. This field is entitled, 

"WATER LEVEL (CM) AFTER T= " 

Z2L The height of the interface above the inital 

position in CM. This field is entitled-, 

"MLD DEPTH (CM) AFTER T= " 

RAD, ANG The top^layer resultant current speed (RAD) in 

CM sec and the direction (ANG) toward -Which 
the current is flowing in degrees true (North is 
0°, East is 90°). This field is entitled, 
"RESULTANT CURRENT SPEED (CM/SEC) . AND DIRECTION . 
(DEG. ) AFTER T= " . . Note that the field' consists 
of numbers of the form RR.DDD. where RR is .'the 
Current Speed and DDD is the Direction. 

RAD2 ,ANG2 The bottom layer resultant current. speed ( RAD 2 ) 

and direction (ANG2) in CM sec arid degrees true. 
The same conventions are sued for RAD2 and ANG2 
as for RAD and ANG. The fiedl is entitled-, "DEEP 
LAY. CURRENT SPEED (CM/SEC) AND DIRECTION (DEG.) 
AFTER T= ". 

S The temperature of the top layer, in degrees Centi- 

grade . The field is labeled "TEMPERATURE OF SUR- 
FACE LAYER (DEGREES) AFTER T= " . It should be 

noted that these 30 by 25 fields are printed in . 
two sections of 15 by 25 each. They are indexed 
properly with row and column indicated so there 
should be no problem with interpretation. 1 
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For DEBUG=1 (debug print desired), the following are 
printed as the first ten columns of the 25 by 30 
array upon initialization: 

The components of the wind field (in msec ^ ) and 
the value of KEY which was used to select the 
option- The first twenty-five lines are. those 
of the XW array (the east-west component of the 
wind). The next twenty-five lines are those of 
the YW array (the north-south component of the 
wind). The print commands are located in sub- 
routine SOI, cards 26,27. 

The components of the initial current field 
(CM sec - ) and the value of KEY which was, used to 

select the option. The first twenty- five,. lines of 
the array are those; of U (the east-west component 
of the top layer of the current) . The -next twenty 
five lines are those of the V array (north-south) 
component of the current in the top layer) . 

Note that the bottom layer is initially at rest. 

The temperature field. in degrees Centigrade. As 
only one option is available at present, no value 
of KEY is printed. 


TABLE A-l 

RANGE OF TYPICAL. VALUES 


PARAMETER 

MIN / 

MAX 

ALPHA 

.995 

. 999 

ALP 

.90 

. 99 

AUS 

5 .-OxlO -4 

5. OxlO -3 

C 

1. 8xl0' 6 

3. 6xl0~ 6 

R 

2. OxlO -3 

4.0x10” 3 

RHOl 

1 . 023 

1.028 h; 

RH02 

1.0 24 

1.029 
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SPECTRO-ANGULAR WAVE MODEL 
FOR THE STUDY OF THE MODIFICATIONS 
OF SURFACE WAVES AT OCEANIC FRONTS 


CONTENTS 
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1 • The Basic Spectro-Angular Wave Analysis/ 

Forecasting Model (The French Model) 


The spectro-angular wave analysis/forecasting model 
used in the present study has been developed in Meteorologie 
Nationale, Paris, France (Gelei and Devillaz, 1969 - 
Devillaz, 1965, a,b). It has been, however, considerably 
modified by the present authors and adapted for the specific 
purpose of this study. This French spectro-angular wave 
model (DSA5) is considered to be one of the .best available 
numerical wave forecasting models. The method of propagation 
and spreading of wave energy makes this model attractive for 
this particular study. 

The basic formula for computation of spectro-angular 
energy density is: 

. Ae 


de 
dt 

with: 


T' 


W, 0-w ,T 


T 


. m 
4 o 


- V (T , 0 ) Ve (T , 0) 


~ ^T^0 


where: 


e = the spectro-angular density 

. ^ 

V(T,0) = the "unit speed" of the component, i.e., 

oT . 

V(T, 0) = , g being gravity 

W = surface wind speed in knots 

0 = a component azimuth 

w = a surface wind azimuth 
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T — the wave period in seconds 
Ty= the spectro-angular increase 
A = the empirical dampening coefficient 
The form of the spectro-angular increase is : 

T*' £t,|0-w|,.wJ = P (T,W, ) S (0-w) 


where : 

P = the spectral increase 

S = the function for angular dispersion, such 'that 


r + 2 


S (0-w) d (0-w) = 1 


JLJ 

2 


From y , the decay coefficient of a component is 
given by: 


1 de 
e dt 


.6 4]xir 4 . .1 

4 4 

g dedu T 



where: 

edTdO = the energy density due to that component 
dedu = the specific water mass 
g = gravity 
A = 18 x 10 -9 CGS 

leaving: 


Am 

o 
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The energy density is computed for 6 different periods 
and 16 directions. The growth and decay,, constants are pre- 
computed into a look-up table in the program. 

2. Model Inputs and Outputs 

The wind speed and direction at each grid point is 
prescribed every 3 to 12 hours, depending upon the avail- 
ability of surface wind analysis/forecasts. In the -present 
program, the wind speed and direction is printed out for 
checking in a code as used in the program (speed in 5-knot 
intervals; direction, in 16 compass points), as well as in 
conventional form. 

A special subroutine in the program permits the intro- 
duction of a cyclone. The radius of maximum winds, the 
maximum wind speed and the direction and speed of movement 
of the cyclone is determined by internal parameters , set 
in the main program. 

As the energy propagation is related to grid size, 
the coefficients AMUL and AMULl in the subroutine CONST. must 
be calculated with the following formula:. 

AMUL = 4.5/grid size in nautical miles. 

Grid. sizes smaller than 30 nautical miles should 
not be used if wind speed is introduced each 3 hours or 
longer time interval, as "holes" in energy, propagation 
will occur in the program. 

Various outputs of wave parameters are provided for in 
the program. First the (average height of the "tenth 


highest wave") can be printed at each grid point. The wave 
period and direction of maximum energy are also printed out in 
field form. In order to display and study the energy spectrum, 
a special subroutine IMPDENS prints the distribution of energy 
by compass directions and by different periods at desired 
output points. 

3 . The Organization of the Numerical Model 

The main program, TYWAVE , reads the. input cards, calls. all 
other subroutines and counts the time. It also stores the 
energy densities from previous time steps. 

The second subroutine, AJUST, is for adjustment (and 
smoothing) the energy spectrum. 

Subroutine ECRI computes and prints the height , 

period and direction of maximum energy and prints out these 
fields. 

The subroutine CALCUL computes the actual spectro- 
angular density fields (R 0) considering decay (using proper 
decay constants from AMOR table) and growth corresponding to 
new wind fields and time (using proper growth constants from 
table CROIS) . 

Subroutine VENTOS calls the wind fields, puts the speed 
and directions into codes used in. the program and prints out 
the wind fields. 

Subroutine CONST contains the growth (CROIS) and decay 
(AMOR) tables and other constants required by the program. 
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Subroutine ADVEC simulates the movement of the "swell” 
(waves from previous time step). The "spreading counters". 
AM, BN and CP are generated by subroutines AVANCX and AVANCY, 
which are called before calling ADVEC. 

Subroutine HURRW creates a moving cyclone. The param- 
eters of the cyclone (grid size, speed and radius of .maximum 
winds and direction and speed of the movement of the cyclone ) 
are given in this subroutine. Subroutine PPAGE , called in 
HURRW subroutine, outputs' the wind fields created by it. 

The constants needed for output of wave spectra are 
created with subroutine SELCPT. The energy density table 
is created in subroutine IMPDENS , which calls a compass 
language subroutine DECRIPC which packs the data into the 
array TABLEAU. The array is printed when the final pass 
through the subroutine IMPDENS is made. 

4 . Symbols Used in the. Program 

Card Input - 

The following parameters are read with the data cards: 


I HIND 

0/1 

0=forecast, l=hindcast (set to 0, 
usually) 

IRO 

0/1 

0=new start, l=use old sea 

state 

IMER3 

0/1 

l=fields read in for every 

3 hours 

IMERb 

0/1 

l=f ields read in for every 

6 hours 

A zero 

in the latter 

: ; two positions means that the 

fields 


will be read in every 12 hours.,*; since the wind fields are inter 
nally generated every three hours, the current convention is 
IMER3=1 and IMER6=Q. 
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ISTOCK 0/1 l=stacking of fields for the 

special points output - not used 
currently 

NSTEP Blank/Value The number of hours of computation 

IDENS3 0/1 Densities for the special points 

are taken every 3 hours , 
if IDENS3 = 1 

IDENS6 0/1 1 = special points' densities 

taken every 6 hours 

Other arrays and variables in the program are : 


AM 

AMOR 

AMUL 

AMULl 

BACK 

BN 

COST 

CP 

CROIS 

EX 

EY 


The longitudinal spreading 

The decay constants 

Coefficient of advection velocity 

A decay coefficient relative 
to grid, used to compute the 
the overall decay 

Background wind speed (not used in 
this model) 

The later spreading counters 
(odd) 

Cosine table for advection 
computations 

The lateral spreading counters 
(even) 

The growth- constants 

Spreading counter in the 
x direction 

Spreading counter in the 
y direction 
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F 

Incrementation value for the 
three spreading counters ; 

AM, BN, CP 

IALT 

A counter, with the value of 
either +1 or -1, which is used 
by the program to consult the 


IDIFM tables 

JAN 

Year 

ICHART 

A counter which triggers output 
of the sea-state graph 

ICOMPM 

Longitudinal spreading .reference 
cycle constant 

ICOMPN 

Lateral spreading reference 
cycle constant (odd) 

IDE 

Same as IDENS 

IDIFM 

■ ■ The" longitudinal ■-spreading ! '‘Gonst ant s 

IDIFN 

The lateral spreading constants 
(odd) 

ID IFF 

The lateral spreading constants 
(even) 

IHED 

. Time coup. ter 

IHER 

. Hour 

INDRFT 

Indrift component of the cyclone 

INTEN 

Maximum wind, speed , knots 

I RANGE (IRAN) 

A counter used to load the fields 
as .a .safety factor 

IRESUL 

The array containing the height of 


the 1/10 highest waves, 
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ISTEP 

IT 

ITETA 

ITETMA 

IVEN 

IVENT 

(IV2) 

JOUR 

MESH 

MOIS 

NR, NW 

P 

Q 

SINT 

RCIRC 

REYE 

RO 

ROCUM 

ROCUM1 


Step counter 

Index counter for the 6 period loop 

Index counter for the 16 direction 
loop 

The array containing the "maximum 
wave direction" 

A counter which causes the next 
set of wind fields to be read 

The two arrays which contain the 
input wind. data, the first being 
wind speed and the second' . the 
wind direction 

Day 

Grid size (nautical miles) 

Month 

Counters 

Grid step constant, used in AVANCX 

Grid step constant, used in AVANCY 

Sine table used for advection 
computations 

Radius of cyclone (max. Winds) 
(nautical miles) 

Radius of the eye 

The spectro-angular density -field 
generated for a given time interval 

The summation of the spectro-angular 
densities for the run in progress 

The summation of the spectro-angular 
densities from the previous run 
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ROPRIM 


ROTECU 


ROTETA 

ROTM 
( IROTM) 

ROTT 

UCNTR ■) 
VCNTER -> 


The summation field of the spectro- 
angular densities used to add the 
advection effects to RO 

The summation of the spectral direc- 
tions, used to find the direction 
of maximum energy 

The "maximum direction" retained 

The maximum period field 

The six period arrays from which 
the ROTM value is chosen 

U and V components for the movement 
of cyclone center (knots) 


VMAX 

XORIGY 
YORIG j 

XBACKA 

YBACK 


Maximum wind speed of the cyclone 
(m/sec) 

X and Y for the original center 
of the cyclone 

X and Y components of the background 
wind speed (knots) 
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WAVE REFRACTION MODEL 


CONTENTS 
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Program 


The Dobson Wave Refraction Model 
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Organization of the Numerical Model 
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listing separately submitted 


The Dobson Wave Refraction Model 


Dobson utilized the theory of Munk and Arthur as the 
basis of his model. Figure 1 gives the definition sketch for 
the wave ray equations. Application of Huygen's principle 
for monochromatic wave propagation yields the following 
relationships: 


Also 


<5c<5 1 


Thus, as t ->• 0 , n -* 6n 0 and 

6a ■' _ Dc 
fit Dn 


~1 Dc 
c Dn 


From this it follows 


n = n.B where B = -rr 
0 b. 


fin = 


6a = 


n 0 6B 


n g SB 
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Using Gauss's integral theories this yields. 

Pet _ 1 dB 

Bn B ds 

which in turn, after applying operators and combinations, 
results in the equation for ray separation as a function of 
time: 

P ,pa „ P 2 c 

Pn l Dt ; " _ 2 

Bn 

After further manipulation, this equation.. can be converted 
to the form 

+ p(t)|| + q(t)B = 0 

dt 

where: ■ 

p(t) = -2 (cos (a) ~ + sin (a) ~) 

and 

2 2 
q(t) = c(sin 2 (a) + 2 sin (a) cos (a) 


+ cos 2 (a) — — j) 

3y^ 

The equation for the curvature of wave rays is then; 

T . 1 . , . ■ 3c ’ , , 3c 

K = e slnW S - cos( “> ay 

The iterative equations which result in the summarized 

solution are as follows: 
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Aa = K g As 


K 

s 

II 

i 

w 

H* 

+ 


1/2 


a i+l 

= a . + 

1 

Aa 




a 

11 

p 

H- 

+ 

a ) 1/2 = 

“i+1 

a . + 

r 

Aa 
■ 2 

X i + 1 

= X. + 
1 

As 

cos (a) 



Y i + 1 

- y . + 

i 

As 

sin (a) 




If further operations and the chain rule are applied, 
the governing equations can be expressed as functions at the 
mean water depth, h(s,y), as follows: 


9c 

9h 

9c 



9X 

9x 

9h 



9 2 c 

3 2 h 

9c 

( ih, 2 

^9x' 

9 2 C 

9x 2 

3x 2 

9h + 

9h 2 


with corresponding expressions for the y. - direction terms. 
The wave celerity has previously been given as: 
c 


= C Q tanh (•—!) 


This form can be readily converted to 


dc 

dh 


ac (1-R c ) 


cR + ah (1-R ) 

c c 


where 


c_ 

x. 
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Applying the chain rule to this equation yields: 

dc 


,2 

d c 


U 


dh“ 


dh 


where 


U 


-2ac R 

c 

[cR + ah ( 1 - R 2 ) J 2 

V 


Substituting into the governing equation, yields the 
final form of these equations: 


u ( H > 2 


+ « (^) 


a 2 c 

dc 

r^h 

ax 2 

dh 

ax 2 

a 2 c _ 

dc 

a 2 h 

3y 2 

dh 

Lay 2 

3 2 c _ 

dc 

a 2 h 

axay 

dh 

axay, 


+ u 


ay 


ah ah 
ax ay 


Dobson uses . the least squares method and a second degree 
equation to describe the depth from the trend surface 


= e. +e-x+e_y.+e.x + e c xy + e,y 

i 2 3 J 4 5 J 6 * 


The deviation. at that grid point is given by: 

Ah = H . . h, 

il ' t 
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His method utilizes twelve points and coefficients chosen to 
minimize the squared deviations . With the least squares 
criterion that: 


; 2 

F (e-^f s 2 1 e 4 / e 6 ^ — 2 ( Ah) 

then there will be a minimum when: 



Dobson uses a matrix approach to this problem and finds 
polynomial solutions from: 


[E] .= [sxy] 


[H] 


where [E] is the coefficient matrix (e^ — - e^) 

[sxy] = [xy] _1 • [uxy] 

[uxy] is a constant, special unit matrix 
[H] is the depth 'matrix 

This yields the derivatives of h with respect to x and y as 
shown: 

3 h , _ • . , 

3x = e 2 + 2e 4 X + e 5^ 


3 2 h 

3x 2 


= 2e . 


3h 

3x3y e 5 


with similar expressions for the y terms, 
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The curvature and wave separation are now representable 
in terms of depth using these coefficients. The equations . 
become 

o (1-R 2 ) 

K = [sin (a) 

cR + ah (1-R ) 

c c 


(e 2 + 2e 4 + e 5 y) -.cos (a) (e 3 + e^x + 2e g y).] 
-2ac (1-R 2 ) 

p(t) .= - — -5- 

cR + ah (1-R ) 

c c 


[cos (a) (e 2 + 2e 4 x + e g y) + sin (a.) (e_ +— e_x + 2e g y) ] 

ac 2 (1-R 2 ) 

q(t ) = 5 __ 

cR + ah (1-R ) 

c c 

{sin 2 (a) [2e 4 + U (e 2 + 2e 4 x + e^y) 2 ] + 

2sin (a) cos(a) [e^ + U ( e 2 + 2e 4 x + e g y) (e 3 +e 5 x+2e g y)] 

+ cos 2 (a ) [2e g + U (e 3 +e 3 x+2e 6 y) 2 ] } 
where x and y are local coordinate values. The p(t) and q(t) 


solutions can then be used to find 3 by: 


(n 
i+1 


[ (t Pi -2)B^[ + 2 (2-t 2 q i )B ; [ n) ] 
(2 + tp^^) 



(N-l ) 
i+1 


(2+t Pi ) 


where -t is the constant increment of time between point's#- ? 
is an error term, and n is the order of approximation. B is 


then used to find the refraction coefficient. 
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To the above equations, the formulae devised by Johnson 
(1948) have been included in order to simulate the refraction 
at wave fronts by currents. His method is based on a version 
of Snell ' s Law: 

C 0 


u + c 
sin (a) sm (3) 


These relationships are displayed in Figure 2. 
Solving for sin (3) we have 


s in ( 3 ) = sin(a)/{l-m sin (a)) 


where m = 


u_ 

c. 


There results, dependent upon the angle of interception., 
changes in the wave length, height and steepness: 


L = 


H 

L 


J 0 


(1-m sin (a) ) 
H, 


M (A 
L 0 


where 


M = [ (cos (a) /cos (3) { (1-m sin (a) ) 6 / (1+m sin (a) ) } ] 
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C 0 = Initial Wave Velocity 

C = Final Wave Velocity Relative 
to Moving Water 


Figure 2. Refraction of Surface Waves by Currents 


(After Johnson) . 



2 . 


Model Inputs and Outputs 


Card 1 
MI, MJ 

IGRCON 

LIMNPT 

NPRINT 

GRID 


DATA CARD SETUP 

(515, 4F10 . 5 , 12 , FIO .,5) 

The maximum values of I and. J 
which define 'the grid. "’Note : I - 

X +1, J = Y + 1 . Maximum 50 by 50 

Grid unit identifier.: 
l=feet 

y-s* ^nautical miles 
3,==meters 

Maximum, number of points to be 
Computed for .each ray. • 

Frequency for. printed output 

•Number of grid units- -per grid 
division 


DCON 

DELTAS 

GRINC 

ITEST 


Card 2, etc 


Conversion factors for depth units 

: --'Minimum length of ; - ^increment -'along 
ray (grid units) 

Length of increment for ray in deep 
water (grid units) 

If. = 0, depth field (Card Two ) is 
. l/iread, If' - .1, .artificial, deep water 
depth field generated 

(16F5.1) 


DEP (I, J) Depth data at grid points. More 

than one card 

XCONST Value of deep "water depth field 

(in units identified by IGRCON) 

Card 3 , (15) 

NOSETS Number of sets of rays of different 

'•^periods . 

Card 4 (18A4) 

Identifying title for each set 

C-,10 


TITL 



Card 5 
NORAYS 
T 

HO 

Card 6 
X, Y 

A 

OUTPUT VARIABLES 
NPT 
WLO 
CO 


(13, IX, 2P6.1) 

Number of rays in each set 
Wave period (seconds) 

Wave height on deep water {feet) 
(3F7.2) 

Coordinates of the starting point 
for each ray 

Initial direction of wave ray, 
measured anti-clockwise from 
positive X direction (degrees) 

(not previously defined) 

The number of the point on the. ray 

Deep water wave length (feet) 

Deep water wave speed 
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LISTING OF INPUT CARDS FOR AN ACTUAL DEPTH GRID;' (Example) 


37 36 1 300 51650,00000 l.OOOOQ .Q05CO ,20000 1 - 

120.011^.0101.0 96.0 9i.o 83.0 ,34.0 30.0 44.0 /2.0 69,0 65,0 65,0 62,0 55,0 5.2.0 
^6,036,0 30,0 15.0 0 -1,0-10,0-10.0-10,0-10,0-10.0-10 ,0-10,0-10 , 0-10 , 0^10-, 0 

-13. 0-10. 0-10. 0-10.0123, 0113, 0110. 0103,0 94.0 37,0 84,0 61.0 74,0 72,0 6 9,0 '6 4 ! c 

63.0 60,0 55.0 54.0 «3,C 36,0 20,0 5.0 0 -1 , 0 -1 0 , 0 -10 , 0 -1 0 , Q-lo , 0 -ToVo -T : 0‘, ; 0 

-10. 0*10. p-lO.C-10. 0-10. 0-10, 0-10. 0-10. 0125, 0115. 0103.0102.0 96,0 9.0,0 37,0 84,0 

81.0 7-7,0 40,0 64.0 60,0 58,0 53,0 49,0 43.5 32,5 6,0 5,0 0 -1 , C -TO , 0-10 . 0 

-10.0 -10, 0-10. C -10. 0-1 0.0 -10. 0-10, 0-10.0-10,0-10. 0-10, 0-10. 0129, 0115. 0101,0 99,0 

96.0 93,0 90,0 73,0 .75,0 73.5 71,0 63.0 5-7.5 56,.-. 54.5-44.5 39 , 0 30.V0 17,0 SSyO 

0 - 1 , 0 - 1 0 , 0 - 1 0 , 0 - 1 0 . 0 - 1 0 , 0 - 1 0 . 0 - 1 0 . C -1 U . 0 - 1 0 . 0 - 1 0 , G - 1 0 , G -1 0 , Q - 1 0 .0 -1 0 , 0 - 10 ,' 0 

l33, 0136,0118, 0108 ,0 93,0 94,0 40,0 85,0 75.0 41,0 63,0 60.0 55 ,Q- 5.1 yO 

35.5 32,0 14. C 6.0 -1.0-10,0-10 .0-10.0-10,0-10.0 -1,0 -0.5 -0,5 -i,o -0,5 0 

-1,0-10,0-10, 0-10. 0135, 0124, 0114, 0109, 0104.0. 96.0 91,0 86.0 74.0 68,0-64,0 59,5 

55.0 49,0 41,5 39.5 36.0 30,0 9,0 0 -1,0-10.0-10,0-10.0-10.0.-1,0 -0 -0 

0 0 G 0 -1,0 -1,0-10, 0-10. 0138, 0125, 0116,0109,0 58,0 94,0 89,0 84,0 

80.0 72,0 64. C 57.0 50.0 44,0 4 0 „q 39,0 34. Q 24,0 lO-.D 0 - 1 0 . 0 -10, 0 - 1 0 , 0 - 1 0 . 0 

-1,0 -0,5 0 .5 l.C 1,5 1.0 0 0 -O.f -1,0 -1.,.0 136 . 0 l3-2 v 0122v,,'0l 0 8 , 0 

102.0 97,0 -94,0.86.0 73.0 66,0 60.0 55.0 49. Q 43, (V 40 . 0 39,0 3 4 , 0 .12 2 -. 0 ^0 : -1 ; C 

-10.0-10,0-10,0 0 0 0 3 ,5 .5 1,0 1,0 i,0 0 -G,'5 -0,5 -G,5 

137.Q130, 0126.0115. niO? , 0 1 G 2 . 0 92,0 76.0 7q.O bS.r 56,0 50.0 44.0 4i, 0 37,5 35 , n 

31.0 l^.O 6,0 -1,0-10,0-10,0 0 0 .5 .5 1,5 2.0 2,0 1,5 1,5 1.0 

0 0 0 -1.0123. 2125,0114. 0109, C1C1-0 96. r. <?i , 0 62.0 6 5 , 0 6 0 . 0 57 , 0 4 9 . 0 

43.0 41,0 37,0 34.5 29,0 l5,C Q -1.0-10.0-10,0 0 l.G 1,0 1,0 Q 0 

1.0 1,5 1.5 1.8 2.0 1.5 0-1.0l32;0122. 0114. 010?. 0104,0 99,0 S9,o 79,0 

72.0 63,0 53,0 46.0 44,3 39,3 35,0 34,0 30.0 5,0 -1,0-10.0-10,0-1,0 0 1 0 

1.0 0 0 0 6,0 1,0 1,5 l.G 2,0 4, 3 0 0 13 6 . 0129 , 0 1 2'0 , Oil 0 , 0 

lC6,o 53,0 64,0 77. C. 69.. 0 .53,0 49,. 0 30.0 1.7,5 17.0 26,0 24,0 13,0 -1,0-10,0-10,0 
-10,0- 0- 0 0 0 0-10,0 0 0 5.0 2,0 1.5 1.0 "",5 1,0 X, 0 

134.0132,0120, 0.11.3. 010,6. O' 96.0. 77,0 60.0 52.0 43,3 31,0 26,5 2o,0 17,5 23,5 19.0 

12.0 -1,0 41,0 0 30 0 0 Q 0 0-10,0 0 7,0 6,0 2,0 

2.0- 1,0 .5 .5133.0132,0123.0113.3105.0 93,0 61,0 4S , 0 35,0 32,0 29,0 26,0 

25.5 24,5 21.0 20.0 10,0-10,0 0 0 30.0 30,0 20,0 0 0 .0 Q 0 

Q 0. 7,0 2.5 1,0 2,0 2,5 1 , 5141 . 0132 .-3 120 , 0U4 , 0 10 6 , G. 62;, 0' 59, O' 36,0 

34.Q 31,0 29.0 2.9,0.27,0 24.0 21,0 13.0 15.0 0 0 31,0 26,5 2 2.5 21,0 0 

0 10,0 5,0 00 0 0 0 0 0 0 0144,0135,0126,0108,0 

95.0 60,0 40,0 36.0 32.3 30,0 33,0 32,0 30,0 29.0 23,0 13,0 20.0 43,0 32,0 33,0 

3l.a 35,0 37.0 42.0 69.0 50,0 47,0 43.0 00 0 0 0 0-10,0 0 

144 , (3136, 0129. 0115.0 96.0 53,0 44 . 0 33.0 33.’0 34.5 34,0 35 , 0 33 , 0 37 , 0 <0 , 0 4 5 , 0 

51.5 52,0 45,0 45.0 45,0-45,0 45,0 45,0 45,0 45,3 <5,0 71,0 71,0 68,0 43,0 34,0 

15.0 0 0 . 0143,0136.0129,0125.0 69,0 46,0 37.0 3i,0 36,5 36,5 37,0 69, G 

40.0 47,0 56.0 .63.0 65,0 45,0 45,0 45.0 45,0 45,0 47,0 53,0 43,0 45,0 43,0 72,0 

83.Q 81,0 67.0 49.0 37,0 26,0 13.0 15,0150.0146.0138,0127,0104,0 59,0 43,0 39,0 
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LISTING OF INPUT CARDS FOR AN ACTUAL DEPTH GRID (Continued) 


37.0 37,0 41. 0 4?. o 47.0 54.0 53.0 55.0 45.0 45.0 36.0 34,0 34.0 36 .(1 30,0 23.0 

31.0 31.3 10.0 29.0 45.0 0-5.0 65.0 60,0 54,0 39.0 35,0 23,0153,0147.0143,0132,0 

107.0 73,0 46,0 41.0 39,0 39,0 40.0 43.0 50,0 54,5.54,0 51.0 5i,0 38.0 32,0 24,0 

19.0 1.9,0 16,0 14.0 13.0 12,0 10.0 7,0 7.r) 3,0 42.0 53. 0 63.0 66.0 93,0 45,0 

t56. 0152, 0144,0124, 0113,0 97.0 52.0 44,5 43,0 42, n <1,0 47,0 49,0 5.2,0 91,0 .51,0 

50.0 3 9 .0 30.0 23.0 20.0 l7,o 14,0 13.0 8,0 8.0 9,0 0 0 : 0 0 -4 8 , 0 

33.0 30,0 40.0 51. 0160,0154, 0147. 0132, Olio, 0 93.0 71.0 52,0 45,0 43,0 44,0 49.0 

49.0 51.0 51.0 53.0 52.0 53.0 29,0 6.0 6,0 6.0 6,0 6.0 6, 0 6,.o 6,0 0 

-1,0-10.0 0 0 42,0 35,0 22,0 21.0168,0159,0150.0139.0122,0107,0 95,0 60.0 

50.0 50,5 50,0 51.0 51,0 51,0 45.0 41.0 46,0 45.5 46,0 0 0 0. 0 0 

0 0 o -1.0-10,0-10,0-10,0 o 0 42,(1 36,0 21,0177* 0-1 6-8 . -0.1- i i-> r 01 -3 7 . 0 

129.0 1 1 4 , 01 03. C 72.0 6-1,0 59.n 55,0 54.5 53,5' 46,0 41,5 41.5 44. 0-'*4 j3, ¥ :0-**5- c 0-;4.2,*: 0 

38.0 36.0 32.0 29.5 24.0 ?C , 0 15,0 9.0 0-10, P-10,0 -1,0 \ U 22,0' 42,0' 47, 0 

tfiO. 0177 , 0162,0147.0137,0123, P 111,0 92. 0 75,0 65 , o 64,0 61.0 57,5 52,5 43,5 47,5 

46.5 46,0 44,0 42.0 39.0 39,0 37,0 34,0 31.5 26, n 19,0 14,0 5,0 0 -1,0 -1,0 

0 ll i 5 23,5 41.0179,0174.0162,0143.0137,0121.0118.0113,0 97.0 85,0 77,0 64,0 

59.5 56,0 4o,(j 50,0 52,0 48.5 47,0 43.0 40.0 39.5 38,5 37,0 34,0 30,0 25,0 13,0 

11.0 0-lO.n -1.0 000 20, 0163. 0176. 0171, 0158, 0143,0133. 0129*0121.0 

tC7 , 0 91,0 ae.O 84,0 71.0 58.0 57.0 57.5 57,5 55.0 51,0 47.0 43,0 41,0 39,5 39,0 

37.0 34,0 29,0 23,5 16,0 5,0 0 -1.0 -1.0 -1.0 Q 0165,0183,0178,3174.0 

163.0149.0133.0107.0110.0104.0 99,0 88.0 78.0 74.0 71.0 65,0 57,0 56,5 56,0 50.5 

46.0 42,0 40,0 39.5 9.0 3«.fl 31.0 27.0 21.0 11.0 3 , 0 - 1 . 0 -1 0 , 0 - 1.0,, 0 -.1 0 ,0 - 1 0 , 0 

700. 0192. 0179. 0168. 0154. 0146.0132.0132. 0122. QUO. 0109,0 92.0 67.0 79,0 73,0 70.0 

65.0 60,0 54.0 52.0 51, C 47 . .0.44,0 43.0 41.0 39,0 65,0 30,0 24,0 17 , o 5,0 0 

-1. 0 -10 . 0 -1 0, 0 -1 0 . 0 20 1, 0 20 0, 0 195, 0 1 84. 0166. 014 4 . pl35 ,0126. 0121, 0.121, 01 16,, 0'i 0 3.0 

93.0 96,0 82,0 79.0 63,0 66.0 59,0 55.0 63.0 50.0 47,0 46.0 44,0' 42, -O 39,0 33,0 

27.0 24,0 11,0 0 -1.0-10,0- 10 , 0-10,0202.020 0,0193,017 9, 0174, 015 6 ,0.1 61, 014 5,0 

132.0125.0119.0108.0 99.0 94.0 39.0 82.0 78.0 68.0 60,0 59, 0 57,0 5i,o 50,0 49.0 

47.Q 45.0 -11,0 39.0 31,0 23,0 17,0 3.0 0 -1.0-10.0-10.0203,0201.0199,0135.0 

179, 0168, 0154. 0150. C132, 0125.0119.0115.0109,0 92. 0 89,0 83.0 76.0 71,0.63,0 60,0 

58.0 54,0 51,0 50.0 49.0 46,0 4l,C 38.0 34,0 23.0 19.0 11,0 0 -1,0-10,0-10,0 

306 .0202, 0198, 0133*. 0173, 0161, 0155, 0144. 0136. 0.133. 0124. 0115, 0110* 0 9 7.0 91,0 6 6,0 

79.0 72,0 63,0 61.0 59.0 54.0 52,0 50.0 49,0 45,0.43,5 42.0 4Q.0 31 , 0 23,0 19,0 

2.0 -1,0-10.0-10.020 9.0203,0200 ,0194.0162,01/3.0164.0150 .01.4.7*01 3.5 ,0127 ,,01.16, 0 

108.0105.0102.0 96.0 65.0 79,0 78,0 73.0 66,0 63. 0 58.5 56,0 52,-5. 5.0,0 -48,0 46,0 

40.0 34,0 .29,0 22.0 12,0 0 -1,0-10, 0212. C2l)4 t pl 96, 0189,0163 ,0.1 7.8 . 01 63,01 60 ,0 

159.0153*, 0149,0126. 0125, 5121. 0111. 0104.0 91,0 82,0 78.0 77.0 74,0 67, 0 62,0 5.9. 0 

56.0 54,0 51,0-48,0 44.0 39,0 30,0 25,0 14,0 2,0 -1.0-10,0216.0208,0204,0194.0 

192.0130, 0172. 0165. 0160 , 0154 .0149, 0131 ,0129.0123, 0115 ,0101 , 0 96,0 'Vi, 0 8-7,0 61,0 

76.0 72,0 65,0 61.0 57,0 55,0 54,0 49,0 47.0 41, 0 33,0 28.0 20,0 5,0 -1,0-10.0 

220 . 0212. 0200 . 0196 . 0192. 0181 .0170 . 0167. 0163.0151. 0147.0134 . 0127.0122, 011.3, Oil?., 0 

104, QicO.O 97.O 87.0 48,0 72,0 71,0 65,0 60,0 56,0 55,0 52,0 47,0 42,0 35,0 30,0 

22.0 20 , Q -1.0-10,0 

3 

COLUMBIA R , PERIOD 8 SECCS, DEEP WATER DIRECTION. E.S.E, 

30.0010.00 

10.0 3,5 110.0 

12.0 .3,5 lio.o 

27.0 3,5 110.0 

COLUMBIA R, PERIOD 10 SECS, DEEP -WATER DIRECTION, E.S.E. 

3 10,00 10,00 

10.0 3.5 90.0 

15.0 3,5 .90,0 

25.0 3,5 90,0 

COLUMBIA R, PERIOD 12 SECS , DEEP: WATER DIRECTION . E,S,E, 

3 12,00 10,00 

10.0 3.5 25.0 

15.0 3,5 25,0 

25.0 3,5 25,0 
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3. 


Organization of the Numerical Model 

The main program is written in FORTRAN IV and is easily 
adaptable to any operating system with a FORTRAN compiler . 

The main program is largely a bookkeeping routine > it 
reads in the depth and wave data, initializes many of the 
necessary parameters and keeps track of which -ray is -being 
computed. it also sets the size of the common block being 
utilized, a description of which follows. 


Variables in Common 


IPT (200) 

Intermediate storage array 

D (12) 

Depth at grid-, points used for 
surf ac e f i ttirig . computations , 
D (K) . (in depth units) 

E (6) 

Coe f f ic ients for the ...equation of 
surface of best fit, E (L) 

Bl, B2 

Values of Beta at the points 
NPT-and (NPT-1) 

CXY 

Wave speed at a point on the 
ray (ft/sec) 

DCDH 

The differential coefficient of 
speed with respect to depth ' (dc/dh) 

DCON 

Conversion .factor .for depth units 

DELTAS 

Minimum length of increment along 
ray (grid units) 

DRC 

Depth at which ‘ref raction com- 
mences (0 . 6 x WLO) 

DTGR : 

Time interval between calculation 
points on the ray 
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DXY 

GRINC 

HO 

IGO, JGO , KGO 
LIMNPT 

NPRINT 

NPT 

PHX, PHY 

RCCO 

RHS 

RK 

SIG 

sk : 

TOP 

WL 

WLO 

SLOP 

DS 

DY2 

DYB1 , DYB2 

ILGO 

V 

WLB 


The calculated depth at the point 
X, Y 

Length of increment for .ray in 
deep water 

Wave height in deep water (feet) 

Branching identifiers 

Maximum number of points to be 
computed for each ray 

Frequency for printed output 

The number of the point on. the ray 

Partial differential coefficients 
of depth with respect to X and Y 
(3h/3x, 3h/3y) 

Ratio of actual wave speed to 
deep water wave speed > (CXY/CO) 

Maximum limit for X on a ray 

Refraction .. coef ficient 

Angular frequency of the wave (a) 

Shoaling coefficient 

Maximum limit for Y on a ray 

Wave length at a point on the ray 
(feet) 

Deep water wave length (feet) 

The slope at frequency of output 

Incremental distance along ray 

The depth at the new. point 

The depths' at the . breaking point . 
and ' reform! ng>p 6 int , respectively 

Intermediate identifier 

Intermediate value in computing 
the curvature 

Wavelength at breaking (feet) 
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SUBROUTINE RAYCON 

This routine controls each individual ray as it 'progresses 
across the grid. Initially it assumes that the wave is still 
in deep water and calculates the second point on the. ray. 

The subroutine DEPTH, is called and the depth at this point is 
found . 

If this depth is greater than the refraction depths-/ then 
the subroutine WRITER is called and the details are printed 
out for that point. If the wave has reached shallow waters 
CURVE is called to calculate the initial value of ray curva- 
ture . 

With the initial point found, the subroutine REFRAG is 
called to find the ...next point on the wave ray. HEIGHT, is 
then called to compute the wave height at that point. If 
the breaking height criterion is not- exceeded , the wave 
details will be printed, depending on the frequency of output, 
or the next point will be computed. 

If breaking is evident, the routines SLOPE and BREAK are 
called. If certain program limitations are met, RAYCON also 
acts to stop the computation of that ray; for instance, if 
the ray has reached the boundary or there is no convergence. 

Local Variables 

ANG The ray angle with respect to the 

X axis (degree's) 

A The ray angle (in radians) 
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SUBROUTINE REFRAC 


This routine solves the refraction equations 

1/2 


Aa 

— 

K s A S 

Ks 

= 

[(Ks) i + (Ks) i+1 

a i+l 

= 

a ^ + Aa 

a 

= 

(a i + a i+l )1/2 

X i + 1 

= 

X i + s cos (o’) 

Y i + 1 

= 

Y i + s sin (a) 


iteratively, to find the next point on the ray. 
Local Variables 


FK 

FKK 

XX, YY 

AA 

DS 

RESMAX 

NCUR 


Curvature at the point NPT. 

Curvature at the point (NPT + 1) . 

Coordinates at the point (NPT + 1) . 

Angle of ray at the point (NPT + 1) 

Incremental distance along .ray (in 
grid points) 

Limiting difference between success 
ive approximations for the new 
curvature . 

Identifier.,' controlling the s.tabili 
of the solution... 
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SUBROUTINE CURVE 


This routine first tests for whether .the wave is in deep 
or shallow water and then computes the local speed using the 
appropriate equation. It then computes the local differential 
coefficients and finds the. curvature at the point by the equation 


K 


a U-R c 2 > 

2 

cR c + crh ( 1 ~R C ) sin(a) [ (e2+2e 4 x+e^y) 


- cos (a) 


(e 3 + e 5 x + 2e g y) ] 


where R =• c/co and e^ through e g are computed fitting coeff- 
icients. 

Local Variables 

Cl Intermediate value of wave speed 

used to solve for the intermediate 
water wave, speed. 

FK Curvature at the point X , Y. 
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SUBROUTINE DEPTH 

This routine determines the local origin for -the -point 
on the ray and then calculates the local coordinates... A test 
is made to see whether the point lies within the same mesh 
square as the previous point. If it does, the new depth is 
calculated using the new local coordinates with. the. coefficients 
computed when the ray first entered the square. Otherwise , 
new coefficients for the surface equation are calculated and 
the depth is found. 

Local Variables 

SXY The special unit matrix for- use 

with a square grid. 

I, J The grid coordinator at the local 

origin for the point at the ray. 
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SUBROUTINE HEIGHT 

This routine calculates the shoaling coefficient and re- 
fraction coefficient using the value of Beta calculated at the 
previous point. These two coefficients are then used to cal- 
culate the local height of the wave . 

A test is made to see if the breaking wave criterion 
has been reached. If so, the subroutine BREAK is called by 
RAYCON which uses the reforming height and wavelength to 
calculate subsequent wave heights. If the wave has not 
broken, it computes the finite difference form of the 
equation of wave intensity to. find the value of Beta at the 
next point. 

Local Variables 

H The wave height at a point on the 

ray. • 

CG The group speed of the wave . 

P, Q Variables used in the calculation 

of Beta. 
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SUBROUTINE SLOPE 

This routine calculates the slope of the surface using 
the depths at each point and then taking a mean slope over the 
frequency output. It thus gives the contour over which the ray 
traverses . 

Local Variables 

SL0P1 The slope of the bottom contour at 

each point . . 

SLOP The- -slope at frequency of output. 


< 
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SUBROUTINE BREAK 


This subroutine. -initialites the wave height to the 
breaking wave height and then calculates the length of the 
breaking zone. It then determines the coordinates arid 'calls 
the routine DEPTH to determine the depth at the reforming point . 

The slope is calculated and from this the reforming wave 
height and length are determined. If the slope is non-reforming, 
the routine will continue to evaluate heights until either a 
reforming slope is encountered or the boundary conditions are 
reached. 


Local Variables 


RUN 

The slope over the breaker zone. 

HA 

The reforming wave height. 

HP 

Wave height prior to reforming. 

LA 

Reforming wavelength. 

WLB 

Length over the breaker zone. 
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SUBROUTINE ERROR 


This routine. estimates a measure of the error by using 
a least squares surface to calculate the depth at the point 
X, Y. It finds the maximum difference between the depths at 
the grid points calculated, using the surface polynomial, and 
the actual data depths, arid then computes the standard devia- 
tion of these differences. The maximum difference is expressed 
as a percent of depth at the point X, Y. 

Local Variables 

DP 1 The computed d e p t h-a t the ... fou r g r i d 

points surrounding the point X, Y. 

DIFMAX The maximum difference between the 

actual depth at a grid point and 
the computed depth. 

FIT The standard deviation of the least 

squares surface. 



SUBROUTINE WRITER 

This routine controls the printed output of the program. 

It also includes various error messages which are printed when 
the error conditions have occured. Also, all intermediate 
stops are handled by this routine. 

Local Variables 

NWRITE Branching identifier to control output 
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SUBROUTINE BAR 


This routine calculates the amount of deformation to 
the wave front crossing an oceanic velocity discontinuity, 
including wavelength and height changes, based on the equations 
described in Section 1. 

Local Variables 


XM 

U 

STP 

XMM 


Defined as u/cxy 
A current velocity 
Resulting steepne s s 
The steepness factor 


Note : 

U is set internally in BAR, slight alterations would 
allow it to be read from a data card. 
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The Formulas Used in the Model 


The heat and moisture exchange program .computes .the 
latent heat exchange .with the following formula (Laevastu, . 
1960) : 

Q = CO. 26 + 0.077V) (0.98e - e ) L. [1] 

q . w •••••" ''t • • • - 

where is latent heat exchange (g cal cm , 24h) , V is 
wind speed (m sec ^) , e y is saturation vapor pressure of the 
sea surface, e is vapor pressure of -the air at- the height 
of 10m, and is latent heat of vaporisation. 

The sensible heat exchange is computed with the following 
formula (Laevastu, 1960) : 


Q, = 39 (0.26 + 0.077V) (T - T ) [2] 

h w a 

-2 

where is sensible heat (g cal cm , 24h) , T^ is sea surface 
temperature (°C) and T is air temperature at the height of 

3 . 


10 m. 

For computation of saturation vapor pressure of the sea 
surface the Gof f-Gratch (1946) formula is used: 


lo ?10 e w 


= '-‘7 .'9 0298 (T /T • ; - 
■ s' wa 

- 1.3816 x id -7 (10 
+ 8.1328 x 10~ 3 (10 


* I) '+ 5-.'02S08 log 1(3 - 

11.344 (1— T /T ) 

wa' s - 

-3.49149 (T /T - 1) 
s' wa 


(T /T ) 
s' wa' 

1) 13 ] 

-1) + log 10 e ws 
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where e is saturation water vapor pressure, e is ..saturation 
w ws 

vapor pressure at steam point temperature/ T is steam point 
temperature (373 . 16°K) and T wa is absolute temperature of 
the sea surface (°K) . 

The sea surface temperature and surface wind are pre- 
scribed (and/or available from other programs). Changes in 
surface air temperature and water vapor pressure r dif€^tniti?al 

O 

are computed with the Amot-Mosby formula with Laevastu-Harding 
(1974) constants. 

3T 

AT - AT o e~ Kt f (f + ,| v r ^)"(i - e Kt ) -[5] 

The same formula is used both for air temperature (T 'V ' and 

cL 

water vapor pressure (e ) computation, whereby AT and T 

cl W 

in the above formula are replaced with Ae and e . The 

w 

symbols used in formula [5] are: AT is T - T , AT ;is 

W a. O 

AT at time Q, t is time (hours) , C and K are empirical con- 
stants (for T : C = 0.10 and K = 0.35 and for e : C = 0.40 and 
a a 

K = 0.50) and 9T w /3r is the rate of change of sea surface 
temperature under the trajectory of the surface wind. 
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2. Model Inputs and Outputs 

In addition to the necessary constants (grid size, -time 
step, etc.), the main inputs required by the model are sea 
surface temperature and surface winds. Synoptic observations 
of surface air temperature and water vapor pressure of the air 
are too inaccurate and too sparce in space and time to use for 
direct computation. Thus, these are computed with the s we 11 
verified Amot-Mosby formula (Laevastu and Harding, 1974). 

The model computes (and prints out) sensible and latent 

-2 

heat exchange in terms of g cal cm in 24 h. Furthermore, 
the resulting changes of sea surface temperature as well as 
the sea surface temperature changes due to advection and eddy 
diffusion can be computed with this model. 
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3. 


Organization of the Numerical Model 

The main program THERMAL reads the constants and other 
data, calls other subroutines, keeps time- and prints out the 
surface temperature field. 

The first subroutine called is SETFLD. It initializes 
the sea surface temperature,, currents, surface wind , and an 
initial guess for the surface air temperature. 

The next subroutine WATER computes the saturation vapor 
pressure of the sea surface and vapor pressure and temperature 
of the surface air. 

The subroutine ANOMALY computes the advection and diffusion 
of sea surface temperature. The final subroutine CALC computes 
the convective transfer of sensible heat .and ..the... transfer 
of latent heat by evaporation . 
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Symbols Used 

in the Program 

A 

Eddy diffusion coefficient 

AA 

Absolute temperature of sea surface 
(°K) 

AK 

TS/AA 

AZA 

Output interval (sec) (data) 

AZZ 

Time counter for output (initialized 0) 

C 

Constant in air temperature computa- 
tion (data) 

DT 

Time step (sec) (data) (also TTT) 

EA 

Vapor pressure of the air (mb) 

EC 

Constant C in vapor pressure compu- 
tation (data) 

EK 

Constant K in vapor pressure compu- 
tation (data) 

EW 

Saturation vapor pressure of sea 
surface (mb) 

EWEA 

(e - e ) 
w a 

EWS 

Saturation pressure of water at steam 
point temperature (data.) 

H 

Time in hours 

K 

Constant K. in air temperature compu- 
tation (data) 

LG 

Grid size (cm) (data) 

NE 

Number of grid points in y direction 

ME 

Number of grid points in x direction 

QE 

-2 

Latent heat exchange (g cal cm , 24h) 

QH 

: v -2 

Sensible heat exchange (g cal cm , 24h) 
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T 


TA 

TE 

TS 

TW 

TWTA 

U 

UU 

UW 

V 

VW 


Time counter (sec) 

Air temperature ( °C) 

End of computation (sec) (data) 

Steam point temperature (data) 

Sea surface tempera. fu re ( ° C ) 

(T - T ) 
w a 

Surface current, u component 
Relative humidity 
Surface wind, :..u component .. 
Surface current, v component 
Surface wind, v component 
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